import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord, BarycentricTrueEcliptic, ICRS
import astropy.units as u
from matplotlib.patches import Circle
import matplotlib as mpl

# 设置支持中文的字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans', 'Arial Unicode MS', 'sans-serif']
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 设置专业天文图表风格
plt.style.use('dark_background')
mpl.rcParams['axes.edgecolor'] = 'white'
mpl.rcParams['axes.labelcolor'] = 'white'
mpl.rcParams['xtick.color'] = 'white'
mpl.rcParams['ytick.color'] = 'white'
mpl.rcParams['grid.color'] = 'gray'
mpl.rcParams['figure.figsize'] = [14, 10]

# 创建黄道上的点 (0°到360°黄经，黄纬0°)
ecliptic_lons = np.linspace(0, 360, 1000) * u.deg
ecliptic_lats = np.zeros_like(ecliptic_lons)

# 将黄道坐标转换为赤道坐标 (J2000历元)
ecliptic_coords = SkyCoord(lon=ecliptic_lons, lat=ecliptic_lats, 
                          frame=BarycentricTrueEcliptic, equinox='J2000')
equatorial_coords = ecliptic_coords.transform_to(ICRS)

# 创建Mollweide投影的天球图
fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection='mollweide')

# 转换为弧度并调整RA范围到[-π, π]
ra_rad = np.radians(equatorial_coords.ra.degree - 180)
dec_rad = np.radians(equatorial_coords.dec.degree)

# 绘制黄道
ax.plot(ra_rad, dec_rad, 'y-', linewidth=2.5, label='黄道')

# 绘制天赤道 (赤纬=0°)
dec_eq = np.zeros(100) * u.deg
ra_eq = np.linspace(0, 360, 100) * u.deg
equator_coords = SkyCoord(ra=ra_eq, dec=dec_eq, frame=ICRS)
ra_eq_rad = np.radians(equator_coords.ra.degree - 180)
dec_eq_rad = np.radians(equator_coords.dec.degree)
ax.plot(ra_eq_rad, dec_eq_rad, 'c--', linewidth=1.5, label='天赤道')

# 标出关键点
# 春分点 (黄经0°, 赤经0h)
vernal_eq = SkyCoord(lon=0*u.deg, lat=0*u.deg, frame=BarycentricTrueEcliptic).transform_to(ICRS)
ax.plot(np.radians(vernal_eq.ra.degree - 180), np.radians(vernal_eq.dec.degree), 
        'ro', markersize=8, label='春分点 (3月21日)')

# 夏至点 (黄经90°, 赤经6h)
summer_sol = SkyCoord(lon=90*u.deg, lat=0*u.deg, frame=BarycentricTrueEcliptic).transform_to(ICRS)
ax.plot(np.radians(summer_sol.ra.degree - 180), np.radians(summer_sol.dec.degree), 
        'ro', markersize=8, markerfacecolor='none', label='夏至点 (6月21日)')

# 秋分点 (黄经180°, 赤经12h)
autumn_eq = SkyCoord(lon=180*u.deg, lat=0*u.deg, frame=BarycentricTrueEcliptic).transform_to(ICRS)
ax.plot(np.radians(autumn_eq.ra.degree - 180), np.radians(autumn_eq.dec.degree), 
        'ro', markersize=8, markerfacecolor='white', label='秋分点 (9月23日)')

# 冬至点 (黄经270°, 赤经18h)
winter_sol = SkyCoord(lon=270*u.deg, lat=0*u.deg, frame=BarycentricTrueEcliptic).transform_to(ICRS)
ax.plot(np.radians(winter_sol.ra.degree - 180), np.radians(winter_sol.dec.degree), 
        'ro', markersize=8, markerfacecolor='none', markeredgecolor='white', 
        label='冬至点 (12月22日)')

# 添加黄道十二宫标记
zodiac = [
    ('白羊座', 0, '♈'),
    ('金牛座', 30, '♉'),
    ('双子座', 60, '♊'),
    ('巨蟹座', 90, '♋'),
    ('狮子座', 120, '♌'),
    ('室女座', 150, '♍'),
    ('天秤座', 180, '♎'),
    ('天蝎座', 210, '♏'),
    ('人马座', 240, '♐'),
    ('摩羯座', 270, '♑'),
    ('宝瓶座', 300, '♒'),
    ('双鱼座', 330, '♓')
]

for name, lon, symbol in zodiac:
    zodiac_point = SkyCoord(lon=lon*u.deg, lat=0*u.deg, 
                           frame=BarycentricTrueEcliptic).transform_to(ICRS)
    x = np.radians(zodiac_point.ra.degree - 180)
    y = np.radians(zodiac_point.dec.degree)
    
    # 确保标签在可见范围内
    if -np.pi <= x <= np.pi and -np.pi/2 <= y <= np.pi/2:
        ax.text(x, y, f'{symbol}\n{name}', 
                fontsize=9, 
                ha='center', 
                va='center',
                color='gold',
                bbox=dict(facecolor='black', alpha=0.5, edgecolor='none', pad=1))

# 添加银河平面作为参考 (简化表示)
galactic_lons = np.linspace(0, 360, 360) * u.deg
galactic_lats = np.zeros_like(galactic_lons)
galactic_coords = SkyCoord(l=galactic_lons, b=galactic_lats, frame='galactic')
galactic_eq = galactic_coords.transform_to(ICRS)
ra_gal_rad = np.radians(galactic_eq.ra.degree - 180)
dec_gal_rad = np.radians(galactic_eq.dec.degree)
ax.plot(ra_gal_rad, dec_gal_rad, 'm-', alpha=0.7, linewidth=1, label='银河平面')

# 设置图表属性
ax.grid(True, linestyle='-', alpha=0.3)
ax.set_title('天球坐标系中的黄道与赤道 (J2000历元)', fontsize=16, color='white')
ax.set_xlabel('赤经 (RA)', fontsize=12)
ax.set_ylabel('赤纬 (Dec)', fontsize=12)

# 设置刻度标签 (转换回度数)
ra_ticks = np.radians(np.arange(-180, 181, 30))
dec_ticks = np.radians(np.arange(-90, 91, 30))
ax.set_xticks(ra_ticks)
ax.set_yticks(dec_ticks)
ax.set_xticklabels([f'{int(x+180)}°' for x in np.degrees(ra_ticks)])
ax.set_yticklabels([f'{int(np.degrees(y))}°' for y in dec_ticks])

# 添加图例
ax.legend(loc='lower right', fontsize=10)

# 添加说明文本
plt.figtext(0.5, 0.01, 
            '黄道是地球绕太阳公转轨道平面在天球上的投影。'
            '黄赤交角约为23.44°，导致四季变化。\n'
            '图中黄色实线为黄道，青色虚线为天赤道。'
            '红色圆点标记了二分二至点，金色符号表示黄道十二宫位置。',
            ha='center', fontsize=11, color='lightyellow', 
            bbox=dict(facecolor='black', alpha=0.7, edgecolor='none'))

plt.tight_layout()
plt.savefig('ecliptic_path.png', dpi=300, bbox_inches='tight', facecolor='#121212')
plt.show()